C
C *** LAST REVISED ON 16-MAR-1988 08:21:37.04
C *** SOURCE FILE: [DL.GRAPHICS.LONGLIB]MAPLIB.FOR
C
C ***********************************************************************
C
C	SOURCE CODE FOR THE WORLD MAP ROUTINES OF THE LONGLIB GRAPHICS
C	LIBRARY.  THESE ROUTINES ACCESS A DATA FILE WITHIN THE LONGLIB
C	GRAPHICS DIRECTORY.  TO ACCESS THIS FILE THE LOGICAL NAME
C	LONGLOC: MUST BE DEFINED TO BE THIS DIRECTORY.
C
C	ROUTINES IN THIS FILE ARE 100 PERCENT FORTRAN-77 COMPATIBLE
C
C ***********************************************************************
C
	SUBROUTINE EARTH3D(R,FLAT)
C
C	WRITTEN BY D. LONG   JULY 1985
C	DRAWS DATA FROM LAND-EDGE MAP OF THE EARTH IN 3d WITH THE OPTION
C	OF A SPHERICAL OR ELLIPSOIDAL EARTH
C
C	R	(R): EARTH RADIUS (IN PLOT UNITS)
C	FLAT	(R): EARTH FLATNESS
C			= 0 FOR SPHERICAL EARTH
C			= 3.3528132E-3 FOR ELLIPSOIDAL EARTH
C
C	READS THE DATA FILE EARTH.DAT IN THE LONGLIB GRAPHICS
C	DIRECTORY.  REQUIRES LONGLOC: TO BE DEFINED AS THIS DIRECTORY.
C	USES FORTRAN FILE UNIT NUMBER 2 FOR READING A FILE.
C
C	USES THE INIT3D 3D PACKAGE BY CALLING PLOT3D
C
	REAL DATIN(10),V(4)
C
C	DEFINE DEGREE TO RADIAN CONVERSION
C
	DEGRAD(X)=X*3.141592654/180.0
C
C	OPEN FILE CONTAINING EARTH LOCATION DATA
C
	OPEN(UNIT=2,FILE='LONGLOC:EARTH.DAT',FORM='FORMATTED',
     $		STATUS='OLD',READONLY)
C
C	BEGIN READING DATA
C
	ID=3
   10	READ(2,1005,END=999) DATIN
 1005	FORMAT(10F8.0)
C
	DO 20 I=1,10,2
		IF(IFIX(DATIN(I)).EQ.0.AND.
     $		   IFIX(DATIN(I+1)).EQ.0)GOTO 20
		IF(DATIN(I+1).EQ.888.8) THEN
C
C	END OF LINE SEGMENT
C
			ID=3
			GO TO 20
		ENDIF
		IF (ABS(DATIN(I)).GT.180.5.OR.
     $		    ABS(DATIN(I+1)).GT.91.) THEN
C	BAD DATA POINT
			IPT=3
			GO TO 20
		ENDIF
C
C	CONVERT INPUT LAT/LONG TO 3D VECTOR
C
		ALONG=DEGRAD(DATIN(I))
		ALAT=DEGRAD(DATIN(I+1))
C
C	COMPUTE EARTH RADIUS ASSUMING ELLIPSOIDAL EARTH
C
		RAD=R*(1.-FLAT*SIN(ALAT))
		ALAT=DEGRAD(DATIN(I+1)-90.)
		CALL SPRECT1(V,ALONG,ALAT,RAD)
		CALL PLOT3D(V(1),V(2),V(3),ID)
		ID=2
   20	CONTINUE
	GO TO 10
  999	CONTINUE
C
	CLOSE(2)
	RETURN
	END
C
C
	SUBROUTINE SPRECT1(V,T,P,R)
C
C	CONVERT FROM SPHERICAL TO RECTANGULAR COORDINATES
C	(LATTITUDE-LONGITUDE TYPE SPHERICAL COORDINATES)
C
C	V	(R):	OUTPUT VECTOR (RECTANGULAR COOR.) V(3)
C	T	(R):	THETA INPUT ANGLE (RADIANS)
C	P	(R):	PHI INPUT ANGLE (CONE ANGLE) (RADIANS)
C	R	(R):	RADIUS INPUT
C
	DIMENSION V(3)
	V(1)=R*COS(T)*SIN(P)
	V(2)=R*SIN(T)*SIN(P)
	V(3)=R*COS(P)
	RETURN
	END
C
C ****************************************************************************
C
	SUBROUTINE LANDMAP(YLAT,XLONG,SLONG)
C
C	WRITTEN BY D. LONG JULY 1985  JPL
C
C	DRAWS A LAND-EDGE MAP USING DATA FROM VECTOR LAND MAP FILE
C	THIS ROUTINE USES A LINEAR PROJECTION, I.E. LAT/LONGS ARE ASSUMED
C	TO BE X,Y POINTS ON A PLANE.
C
C	YLAT	REAL	SCALE FACTOR OF LONGITUDE DEGREES/INCH
C	XLONG 	REAL	SCALE FACTOR OF LATITUDE DEGREES/INCH
C	SLONG	REAL	LEFT-MOST LONGITUDE (-180 TO +180)
C
C	ASSUMES THAT LONGLOC HAS BEEN ASSIGNED THE DIRECTORY OF THE
C	EARTH.DAT FILE STORING THE EARTH LAND-EDGE MAP.  USES FORTRAN
C	FILE UNIT NUMBER 2.
C
C	NOTE: SINCE MAP DATA HAS ITS ORIGIN FROM -180 TO 180, IF WE WANT
C	THE MAP TO START FROM SOME OTHER POINT WE NEED TO DO SOME FANCY
C	TRICKS
C
	REAL DATIN(10)
	LOGICAL LAST,PENUP
C
C	OPEN DATA FILE
C
	OPEN(UNIT=2,FILE='LONGLOC:EARTH.DAT',FORM='FORMATTED',
     $		STATUS='OLD',READONLY)
	IPT=3
	LAST=.FALSE.
	X=2.*XLONG
C
C	READ THE DATA
C
   10	READ(2,1005,END=999) DATIN
 1005	FORMAT(10F8.0)
C
C	CHECK THE STORED DATA COMMAND
C
	DO 20 I=1,10,2
		IF (IFIX(DATIN(I)).EQ.0.AND.
     $ 		    IFIX(DATIN(I+1)).EQ.0) GOTO 20
		IF (DATIN(I+1).EQ.888.8) THEN
C
C	 END OF LINE SEGMENT
C
			IPT=3
			CALL PLOT(X,Y,IPT)
			PENUP=.TRUE.
			LAST=.FALSE.
			X=XLONG*2.
			GO TO 20
		ENDIF
		IF (ABS(DATIN(I)).GT.180.5.OR.
     $		    ABS(DATIN(I+1)).GT.91.) THEN
C	BAD POINT
			IPT=3
			CALL PLOT(X,Y,IPT)
			PENUP=.TRUE.
			LAST=.FALSE.
			X=XLONG*2.
			GO TO 20
		ENDIF
		PENUP=.FALSE.
		IF (DATIN(I).LT.SLONG) THEN
C
C	LEFT OF DESIRED LONGITUDE
C
			IF (LAST) THEN
C
C	CLIP INTERSECTION LINE
C
			    X1=(DATIN(I)+360.-SLONG)/360.*XLONG
				Y1=(DATIN(I+1)+90.)/180.*YLAT
				Y2=Y1-(Y-Y1)*X1/(X-X1)
				CALL PLOT(X,Y2,IPT)
				CALL PLOT(X1,Y2,3)
				CALL PLOT(X1,Y1,2)
				X=X1
				Y=Y1
			ELSE
			    X=(DATIN(I)+360.-SLONG)/360.*XLONG
				Y=(DATIN(I+1)+90.)/180.*YLAT
				CALL PLOT(X,Y,IPT)
			ENDIF
			IPT=2
			LAST=.FALSE.
		ELSE
C
C	RIGHT OF DESIRED LONGITUDE
C
			IF (LAST) THEN
				X=(DATIN(I)-SLONG)/360.*XLONG
				Y=(DATIN(I+1)+90.)/180.*YLAT
				CALL PLOT(X,Y,IPT)
			ELSE
C
C	CLIP INTERSECTION LINE
C
				X1=(DATIN(I)-SLONG)/360.*XLONG
				Y1=(DATIN(I+1)+90.)/180.*YLAT
				Y2=Y-(Y1-Y)*X/(X1-X)
				CALL PLOT(X,Y2,IPT)
				CALL PLOT(X1,Y2,3)
				IPT=3
				IF (PENUP) IPT=3
				CALL PLOT(X1,Y1,IPT)
				X=X1
				Y=Y1
			ENDIF
			IPT=2
			LAST=.TRUE.
		ENDIF
   20	CONTINUE
	GO TO 10
  999	CONTINUE
	CLOSE(2)
C
	RETURN
	END
C
C
C ****************************************************************************
C
	SUBROUTINE POLARMAP(X0,Y0,RAD,ANG)
C
C	WRITTEN BY D. LONG JULY 1985  JPL
C
C	DRAWS A LAND-EDGE MAP USING DATA FROM VECTOR LAND MAP FILE
C	THIS ROUTINE USES A POLAR PROJECTION, I.E. LATTITUDE IS ASSUMED
C	TO BE A LINEAR RADIUS FROM THE POLE.
C
C	X0,Y0	REAL	ORIGIN OF THE POLE
C	RAD	REAL	RADIUS OF EQUATOR
C			> 0 NORTHERN HEMISPHERE
C			< 0 SOUTHERN HEMISPHERE
C	ANG	REAL	ANGLE OF PRIME MERIDAN FROM HORIZONTAL (DEG CCW)
C
C	ASSUMES THAT LONGLOC HAS BEEN ASSIGNED THE DIRECTORY OF THE
C	EARTH.DAT FILE STORING THE EARTH LAND-EDGE MAP.  USES FORTRAN
C	FILE UNIT NUMBER 2.
C
	REAL DATIN(10)
	LOGICAL LAST
C
C	OPEN DATA FILE
C
	OPEN(UNIT=2,FILE='LONGLOC:EARTH.DAT',FORM='FORMATTED',
     $		STATUS='OLD',READONLY)
	DTR=3.1415926/180.0
	AR=ABS(RAD)/90.0
	IPT=3
	LAST=.FALSE.
C
C	READ THE DATA
C
   10	READ(2,1005,END=999) DATIN
 1005	FORMAT(10F8.0)
C
C	CHECK THE STORED DATA COMMAND
C
	DO 20 I=1,10,2
		IF (IFIX(DATIN(I)).EQ.0.AND.
     $		    IFIX(DATIN(I+1)).EQ.0) GOTO 20
		IF (DATIN(I+1).EQ.888.8) THEN
			IPT=3
			GO TO 20
		ENDIF
		IF (ABS(DATIN(I)).GT.180.5.OR.
     $		    ABS(DATIN(I+1)).GT.91.) THEN
			IPT=3
			GO TO 20
		ENDIF
		IF (SIGN(1.,DATIN(I+1)).NE.SIGN(1.,RAD)) THEN
			IPT=3
		ELSE
			ALONG=DATIN(I)
			ALAT=DATIN(I+1)
			R=(90.0-ABS(ALAT))*AR
			A=DTR*(ALONG*SIGN(1.,RAD)+ANG)
			X=COS(A)*R+X0
			Y=SIN(A)*R+Y0
			CALL PLOT(X,Y,IPT)
			IPT=2
		ENDIF
   20	CONTINUE
	GOTO 10
  999	CONTINUE
	CLOSE(2)
C
	RETURN
	END
C
C *********************************************************************
C
C THE FOLLOWING ARE ROUTINES USING THE CIA LAND/SEA BIT MAP FILE
C SEE THE COMMENTS IN THE BITMAP AND LNDSEA ROUTINES
C
C *********************************************************************
C
	SUBROUTINE BITMAP(CLONG,XLONG,YLAT,NXPIX,NYPIX,ILAND)
C
C	BITMAP USES THE FILE "LNDSEA1.DAT" TO DRAW AN EARTH LAND MASK WHERE
C	THE LAND SURFACE OF THE EARTH IS DRAWN USING HORIZONTAL LINES
C	ON A RECTANGULAR PROJECTION OF THE EARTH'S SURFACE (SEE LANDMAP).
C	A GRID OF PIXELS IS USED TO LOCATE POINTS WHICH ARE "LOOKED UP"
C	IN LNDSEA1.DAT (WHICH IS INTERPOLATED) TO DETERMINE IF EACH PIXEL
C	IS LAND/SEA.  USING HORIZONTAL SCANS, A PLOT LINE IS GENERATED OF
C	THE LAND/SEA AREA WHICH IS OUTPUT WHEN THE TYPE CHANGES.  NOTE
C	FOR HIGH DENSITY PLOTS THIS OPERATION IS VERY SLOW DUE TO THE HIGH
C	NUMBER OF FILE READS OF LNDSEA1.DAT REQUIRED.
C
C	INPUTS:
C		CLONG	(R)	LEFT SIDE LONGITUDE OF MAP (DEGREES -180,180)
C		XLONG	(R)	LENGTH OF LONGITUDE AXIS (INCHES)
C		YLAT	(R)	LENGTH OF LATTITUDE AXIS (INCHES)
C		NXPIX	(I)	NUMBER OF PIXELS IN X DIRECTION
C		NYPIX	(I)	NUMBER OF PIXELS IN Y DIRECTION
C		ILAND	(I)	LAND/SEA FLAG
C				 = 0 PLOT LAND
C				 = 1 PLOT SEA
C	OUPUTS: (NONE)
C	CALLS:
C		LNDSEA
C		PLOT
C
	LOGICAL PENUP,LAND,OPTION
C
	XPIX=NXPIX
	YPIX=NYPIX
	CCC=XLONG*(CLONG+180.)/360.
	OPTION=.TRUE.
	IF (ILAND.NE.0) OPTION=.FALSE.
C
	DO 100 ILINE=1,NYPIX
		LINE=NYPIX-ILINE+1
		RLAT=(90.-FLOAT(LINE-1)/YPIX*180.)
		PENUP=.TRUE.
		Y=(YPIX-FLOAT(LINE-1))*YLAT/YPIX
		X0=0.
		Y0=Y
		DO 50 I=1,NXPIX
			RLON=(FLOAT(I-1)/XPIX)*360.-180.	
			X=(FLOAT(I-1)*XLONG/XPIX)
			IF(RLON.LT.0.) RLON=RLON+360.
			IMAP=LNDSEA(RLAT,RLON)
			LAND=.TRUE.
			IF (IMAP.EQ.0) LAND=.FALSE.
			IF (((.NOT.LAND).AND.OPTION).OR.
     $			    ((.NOT.OPTION).AND.LAND)) THEN
				IF (PENUP) THEN
					X0=X
					Y0=Y
					PENUP=.FALSE.
				ENDIF
			ELSE
				IF (.NOT.PENUP) THEN
					X1=X0-CCC
					Y1=Y0
					X2=X-CCC
					Y2=Y
					IF (X1.LE.0.) THEN
					    IF (X2.LT.0.) THEN
						X1=X1+XLONG
						X2=X2+XLONG
						CALL PLOT(X1,Y1,3)
						CALL PLOT(X2,Y2,2)
					    ELSE
						X1=X1+XLONG
						CALL PLOT(0.,Y1,3)
						CALL PLOT(X2,Y2,2)
						CALL PLOT(X1,Y1,3)
						CALL PLOT(XLONG,Y2,2)
					    ENDIF
					ELSE
					    IF (X2.GT.0.) THEN
						CALL PLOT(X1,Y1,3)
						CALL PLOT(X2,Y2,2)
					    ELSE
						X2=X2+XLONG
						CALL PLOT(0.,Y2,3)
						CALL PLOT(X1,Y1,2)
						CALL PLOT(X2,Y2,3)
						CALL PLOT(XLONG,Y1,2)
					    ENDIF
					ENDIF
					PENUP=.TRUE.
				ENDIF
			ENDIF
   50		CONTINUE
		IF (.NOT.PENUP) THEN
			CALL PLOT(X0,Y0,3)
			CALL PLOT(X,Y,2)
			PENUP=.TRUE.
		ENDIF
  100	CONTINUE
	RETURN
	END
C
C ************************************************************************
C
	INTEGER FUNCTION LNDSEA(ALAT,ALON)
C
C	THIS ROUTINE RETURNS A FLAG INDICATING IF THE POINT
C	(ALAT,ALON) IS EITHER ALL LAND OR ALL SEA.
C
C	INPUTS:
C		ALAT	(R)	LATTITUDE (DEGREES)
C		ALON	(R)	LONGITUDE (DEGREES)
C	OUTPUTS:
C		LNDSEA	(I)	LAND/SEA FLAG
C				  = 0 FOR OCEAN
C				 <> 0 FOR LAND
C	CALLS: (NONE)
C
C	NOTE: USES FORTRAN FILE UNIT #1
C
 	INTEGER*4 ALLSAM(481),CIADAT(4,120),BIT(32)
	EQUIVALENCE (ALLSAM(2),CIADAT(1,1))
	DATA IOPEN/0/,LU/1/,IRECLA/0/
C
C	OPEN BITMAP IMAGE FILE OF LAND/SEA
C
	IF (IOPEN.EQ.0) THEN
		IOPEN=1
	        OPEN(UNIT=LU,FILE='LONGLOC:LNDSEA1.DAT',
     &		 ACCESS='DIRECT',STATUS='OLD',READONLY,
     &		 MAXREC=648,ERR=99)
C
C	LOAD DECODING BITMAP
C
		DO 10 I=1,32
			IZ=0
			BIT(I)=IBSET(IZ,32-I)
10		CONTINUE
	ENDIF
C
C   ENSURE VALID INPUT ANGLES
C
	BLON=AMOD(ALON,360.)
	IF (BLON.LT.0.0) BLON=BLON+360.
	BLAT=AMOD(ALAT+90.0,360.)
	IF (BLAT.GT.180.0) BLAT=BLAT-180.0
	IF (BLAT.LT.0.0) BLAT=BLAT+180.0
C
C   SET UP FILE POINTERS
C
	JLAT0=IFIX(BLAT)/10
	IF (JLAT0.EQ.18) JLAT0=17
	JLON0=IFIX(BLON)/10
	IF (JLON0.EQ.360) JLON0=1
	IREC=((JLAT0)*36)+JLON0+1
C
C   SEE IF A NEW RECORD MUST BE FETCHED
C
	IF (IREC.NE.IRECLA) THEN
		READ(LU,REC=IREC,ERR=99) ALLSAM
		IRECLA=IREC
	ENDIF
C
C   CHECK FOR 'ALL LAND' OR 'ALL SEA'
C
	IF (ALLSAM(1).NE.-1) THEN
 		LNDSEA=ALLSAM(1)
		RETURN
	ENDIF
C
C	COMPUTE ADDRESS IN BUFFER
C
	JLAT=INT((BLAT-FLOAT(10*JLAT0))*12)+1
	IF (JLAT.EQ.121) JLAT=120
	JLON=INT((BLON-FLOAT(10*JLON0))*12)+1
	IWORD=(JLON+29)/30
	IBIT=JLON-(30*(IWORD-1))
C
C	EXTRACT BIT
C
	LNDSEA=IAND(CIADAT(IWORD,JLAT),BIT(IBIT))
	IF (LNDSEA.NE.0) LNDSEA=1
	RETURN
C
C	FILE ERRORS
C
99	WRITE(*,45)
45	FORMAT(' *** ERROR OPENING/READING LNDSEA1.DAT FILE ***')
	RETURN
	END
